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Abstract 


In this paper, we introduce second derivative multistep collocation meth- 
ods for the numerical integration of ordinary differential equations (ODEs). 
These methods combine the concepts of both multistep methods and col- 
location methods, using second derivative of the solution in the collocation 
points, to achieve an accurate and efficient solution with strong stability 
properties, that is, A-stability for ODEs. Using the second-order deriva- 
tives leads to high order of convergency in the proposed methods. These 
methods approximate the ODE solution by using the numerical solution in 
some points in the r previous steps and by matching the function values 
and its derivatives at a set of collocation methods. Also, these methods 
utilize information from the second derivative of the solution in the colloca- 
tion methods. We present the construction of the technique and discuss the 


analysis of the order of accuracy and linear stability properties. Finally, 
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some numerical results are provided to confirm the theoretical expecta- 
tions. A stiff system of ODEs, the Robertson chemical kinetics problem, 
and the two-body Pleiades problem are the case studies for comparing the 


efficiency of the proposed methods with existing methods. 


AMS subject classifications (2020): 65R20. 


Keywords: Collocation; Linear stability; Ordinary differential equation; 


Second derivative methods. 


1 Introduction 


This paper is devoted to introducing a new class of collocation methods 
by using the second derivative of the solution for the numerical solution of 


ordinary differential equations (ODEs) in the form 


y' = f(y), t€ [to, 71, 


y(to) = Yo, 


(1) 


where the function f : R¢ —> R¢ is sufficiently smooth, yo € R@ is a known 
initial value, and the dimensionality of the system is shown by d. 

One-step collocation methods for ODEs have been analyzed by many au- 
thors (see [3] and references therein contained). As it is well known, the 
collocation methods is a technique that the given equation is transformed to 
an algebraic equation by representing the solution as a combination of basic 
functions belonging to a chosen finite dimensional space, usually a piece- 
wise algebraic polynomial, which enforces the integral equation at the chosen 
collocation points. 

Multistep methods play a crucial role in the numerical integration of 
ODEs due to their ability to obtain the approximate solution of a high or- 
der of convergence. Lie and Norsett [14] first introduced the idea of multi- 
step collocation methods. In comparison with classical collocation methods, 
these methods depend on more parameters while computational cost does 
not exceed them, and these methods have high convergence order and strong 


stability properties. 
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In comparison with one-step methods, in view of computational cost and 


precision of approximation, the multistep methods are more efficient [13]. 


Using higher order derivative can be a helpful way to achieve the inte- 
gration methods will be more stable with high order convergency in both 
onestep and multistep methods. Firstly, a class of second-order derivative 
formulas is created, and the stability of these formulas was investigated by 
Enright [7]. Following this work, successful implementation of second-order 
derivatives has been developed in many literatures [4, 5, 6]. Also, second 
derivative two-step collocation methods have been recently constructed in 
[8]. These methods and all other traditional second derivative methods can 
be considered as special cases of the class of second derivative general linear 


methods (see, for instance, [1, 2, 4]). 


We aim to provide a second derivative extension of the multistep colloca- 
tion methods (SDMCMs) that use the second derivatives of the solution in 
the collocation points. These methods combine the advantages of both mul- 
tistep collocation methods and second derivative methods to obtain methods 
with higher order accuracy and desired stability properties. In other word, 
we derive a special case of multistep collocation methods, which fixed num- 
ber r of previous time steps and first and second derivative of the solution 
in m collocation points are used in the approximation of the solution in ev- 
ery subinterval. The advantages of these methods are their ability to handle 
stiff ODEs and their superior accuracy compared to multistep collocation 


methods while the computational cost does not exceed. 


Theoretical concepts of new proposed methods, SDMCMs, is discussed 
in details in Section 2. In Section 3, the order condition of the proposed 
method is obtained. Section 4 is dedicated to analyzing the linear stability 
properties of SDMCMs with respect to the linear basic test equation. The 
proposed methods with two and three steps and one and two collocation 
parameters are constructed in Section 5. Finally, Section 6 presents several 
numerical examples in order to validate the effectiveness of proposed methods 


and experimentally verify the order of convergence. 
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2 Construction of method 


Let us divide the interval I := [to, T] to the finite set of subintervals [tn, tn+1] 


for n = 0,1,...,N —1 by introducing a uniform mesh h in the form 
Tn, = {to + nh| n=0,1,...,N, h>0, Nh=T — to}. 


Also, consider the abscissa vector c = [c1,€2,---,Cm]” and define the colloca- 
tion parameters by ty, ; := tn +h. The multi-step continuous approximation 
to the solution of (1) is defined by 


Plin + 8h) =~ 9e(s)om-e +h (vy(8)F(Pltns))) 
k=0 j=l 


Yn+1 =P(tn41), 


with s € (0,1] and the function g is the second derivative of unknown func- 
tion y(t) and is defined by g(-) = f’(-)f(-). Also, for the implementation 
of the method, the approximate solution in the initial interval [to,t-—1] is 
needed where can be computed by an arbitrary method with sufficiently 
high order convergency . The method is defined by the coefficient functions 
ye(s), vy(s), x; (s), 7 =1,...,m, k=0,1,...,—1, which are polynomials 
of degree 2m +r-—1. The interpolation conditions for P(t, + sh) at the 


points t,_;, lead to 
P(tn—i) = Yn-1; 4=0,1,...,r—1, (3) 


and collocating the both sides of (1) and its derivative in the collocation 


points ¢,,; is written in the form 


P" (ina) = TP tia); Ping) 2 o( Pas); 1=1,2,...,m. (4) 


The coefficient polynomials of method are determined by conditions (3) and 
(4), which lead to 


vr(—t) = bx, H;(-7) =0, x;(-2) = 0, (5) 
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and 
yy (ci) =0, bj (ci) = 43, xG(ed) = 0, (6) 
vylci) =0, OF (ci) =90, x7 (ci) = 4:5, 

where 6;;, 7,7 = 1,2,...,m, is the usual Kronecker delta. The construction 


of these polynomials is obtained by the Hermite—Birkhoff interpolation [18]. 
To find the unique polynomial that satisfies these conditions, the polynomials 


pr(s), W;(s), and x;(s) can be formulated in the form 
2m+r-1 si 2m+r—-1 si 2m+r—-1 st 
r ; 
pes)= S) OMS, wls)= YO wa, xal)= DO xP. 
i=0 , i=0 , i=0 , 
(7) 
Now, the result of setting s = —k, k = 0,1,...,r —1, in (7) and setting 
8$=c, t=1,2,...,m in the first and second derivatives of the polynomials, 


leads to a linear system for coefficient of the polynomials. The coefficient 


matrix A € RE™+")x@m+tr) is given by 


1 0 0 O avs 0 
, @P ey. Ger (ayers 
1! 2! 3! “tt (Qm+r—1)! 
1 ( r+1)! ( r+1)? ( r+1)3 (erly rt 
1! 2! 3! “8 (Qm+r—1)! 
ct 3 epene? 
A= 0 1 1 2r “ot (Qm+r—2)! 
a Pe Spee 
0 1 To Grima 
ct eamtr—s 
0 0 1 Tt (mtr—syl 
a pres 
0 0 1 7 (amir—3)l 


and the vectors in the right hand of linear system of equations are defined 
by ul eR, k=1,2,...,r, vl ER™, fj =1,2,...,mas 
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By solving the systems 


All = ful**4 0,,,0n)7,&=0,1,...,r-1, 
Av! =[0,,v!0nJ7, 7 =1,2,...,m, (8) 
Aylil — (0,-,Om, ve)", J = 1,2,...,m, 


the coefficients of the polynomials are determined. 


Now, the conditions for applying the poising condition [9] are established. 
Poising condition is a criterion that determines whether the given set of data 
points can be interpolated using the Hermite—Birkhoff interpolation. By us- 
ing this condition, it can be shown that the solution of the interpolation prob- 
lems can be determined uniquely. For this purpose, the following theorems, 
which are needed to check the uniqueness of the solution of the Hermite— 


Birkhoff interpolation, are mentioned [9]: 


Let & and n be natural numbers, and define the matrix 
B= |le,;||, @=1,...,4, 7 =0,1,...,n—1, 
as a matrix with k rows and n columns having elements 
€;4j =0 or 1, 


which are such that 
S- E55 =n. 
tj 


We shall also assume that the matrix F has no zero rows. Also, suppose that 


{x;}*_, is an increasing real number as 
Ly <%Q2< +++ Xp. 
Also, the set of ordered pairs is considered by 


e={(i,j) | ej = 1. 
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The interpolation problem is described by the reals x; and the “incidence 


matrix” F in the form 
f(x) =y for (4,7) Ee. (9) 


The matrix FE has a structure where each row corresponds to a different in- 
terpolation point, and the columns represent the derivatives at those points. 
It is appropriate to see (9) from the point of view of a Hermite—Birkhoff inter- 


polation problem which we show it in an abbreviated form by HB-problem. 
Definition 1. When the conditions 
P(t) €tm-1, P(«x;)=0, for all (i,j) € E, 


are established, then the HB-problem (9) is poised, and as a result, P(x) = 0, 


equivalently, the matrix E is called poised if the related interpolation problem 
(3) 


has a unique solution for any set of constants y;"’, in which the problem is 


independent of choosing of the ordered interpolation points x, 2%2,..., 2x. 
Define 
k 1 
m=) 6 M=> iy, 91=0)1,n-1 
i=1 j=0 


It is shown by Schoenberg [18] that a necessary condition for poising of F is 
M,>1+1, 1=0,1,...,n—1, 
and these inequalities are recognized as Polya conditions. 


Definition 2. Let the incidence matrix EF be with k rows. Let f; be the 
column index of the first one that appears in the row 7. Moreover, E is called 
a pyramid matrix if, for each i, €;; = 1 implies ¢,;, = 1 for fi < j’ < j, 
and there is some value of 1 < i < k so that f; > fo >--- > fi and 


fi < figa S++ < fee 


The necessary condition for poising the matrix E with respect to the 
ordering points 71 < @2 < ++: < xm, is declared by Ferguson [9] in the next 


theorem. 
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Theorem 1. If EF is a pyramid matrix with k rows, satisfying the Polya 


conditions, then F is poised with respect to the ordering 71 < xg < ++: < &x. 


Now, we show that our interpolation problems (5)—(6) have unique so- 
lution or equivalently the conditions hold true for the given interpolation 


problem. For these problems, the interpolation points can be considered by 
—r+1<-r4+2<-:-<-1<0<c) <cg < +--+ < cm. 


So, the matrix £ with m+r rows and 2m+r columns can be defined by 


[1000---0] 
1000---0 

E=1]1000 
0110 
0110 0 

Thus, we have 
Mo=Tr, M=mM, M2=mMm, ms =), ,Mn-1 = 0, 
Mo =r, M,=r+m, M3 =2m+r, ..., My, =2m+r=n. 


By considering the numbers M;, since 
Mo >=1, M,>=2, M3 >=3,..., Mn_-1 >=n. 
the Polya condition is satisfied for the matrix E. Also, by considering 


fi =9, fo = 9, aay fr = 9, froi = 1, Beas 9 frim =1, 


according Definition 2, one can see that the matrix E’ is a pyramid matrix. 
Thus, by Theorem 1, it is concluded that FE is poised with respect to the 
ordering points —r+1<—-r+2<---<-1L<0<aq <@<--:<cp. Thus, 


satisfying these conditions guarantees a unique and smooth interpolating 
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polynomial that satisfies the given data (5)-(6). In the other word, the 


matrix A is nonsingular. 


When a nonlinear ODE is integrated by implicit (first derivative) meth- 
ods, the approximation of the stage values leads to solving a nonlinear alge- 
braic system of equations, which can be solved by Newton’s iterative methods. 
Thus, the Jacobian matrix Of /Oy is usually computed. Therefore, computing 
the second derivative function g := (Of /Oy)f(y) does not impose any addi- 
tional computational cost. It should be mentioned that in the application 
of these methods, the Jacobian matrix 0g/Oy is approximated by (Of /Oy)?, 


which is a piecewise constant approximation. 


3 Continuous order conditions 


In this section, we investigate continuous order conditions for the method 
(2). We know that the collocation polynomial P(t, + sh) provides a uniform 
approximation of order p to y(t, + sh), s € [0,1]. For this end, the approx- 
imate solution P(t, + sh) is replaced by y(t, + sh) in (2), where y(t) is the 
exact solution of (1). Then by subtracting the both sides of the obtained 


relation, the discretization error is obtained in the form 


€(tn + sh) =y(tn + sh) — y yr(s)y(tn—K) — hye 05(8) f(y(tn.3)) 
. = j=l (10) 
tg > X;(8)9(y(tn3)): 


where s € [0,1] and n = 1,2,...,N —1. By expanding local discretization 
error in Taylor series around the point t, and collecting terms with the same 


powers of h we have the following theorem. 


Theorem 2. Assume that the function f(y) is sufficiently smooth. Then 


the method (2) has uniform order p if the following conditions are satisfied: 
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keels) +> v(s) = 5, (11) 


oo onl) + (+09 ag + xsl) i) — 


where s € [0,1], i = 2,3,...,p. Moreover, the local discretization error (10) 


takes the form 


e(tn + sh) = h?*1C,(s)yt) (tn) + O(APt?), (12) 
in which , 
gPtl c (—k)Pt} 
a +i) d. +i” @) 
(13) 


The set of order conditions (11) leads to a linear system of p + 1 equa- 
tions in 2m +r unknowns. The unknowns of this system are the coefficient 
polynomials of the proposed method. 

Therefore, considering p at most equal to 2m +r—1 guarantees that the 


linear system (11) is compatible. This leads to the following result. 


Remark 1. The maximum attainable uniform order of convergence for the 


new method (2) is 2m+r-—1. 


In the following, the necessary conditions for zero-stability of the new 
method are investigated. For this end, the method is applied to the simple 


ODE 7’ = 0, where the recurrence relation is obtained in the form 


r-1 
Yn4+1 = S On Yn—k) 
k=0 


where 0; = yx(1). Thus, the characteristic polynomial of this recurrence 


relation is 
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r-1 
k 
Pw) = >5 nw, 
k=0 


Therefore, the constructed method is zero-stable if and only if -—1 < 6; < 1 
fori =1,2,...,r—1. 


4 Linear stability analysis 


We now focus our attention on the linear stability properties of the SDMCMs 
(2) with respect to the standard basic test equation of Dahlquist 


y’ = AY, (14) 


where A is a complex number with negative real part. Let us consider the 


vectors and matrices 


A= [dj (cM ER, A= [xi(c) 1 ER, 
vu? = [yi(1)]? € R™, w? = [yi(1)]7 € R™, 


(15) 
(P()i2 = leilelER™,  f=0,1,...,.r-1, 1=1,2,...,m, 


P(tn + ch) =[P(tn + ch) [Phar 


Theorem 3. Applying the new method (2) to the basic test equation (14) 


leads to the recurrence relation in the form 


Yn+1 Mi (z) 0 Yn 
= ; (16) 
y? I, 0,1 yi, 


where z := Ah and 
Mi1(z) = g(1) + (zv7 + 22w7) Q-y(c) € R'*", 
Q=Im—zA—22Ae€R™™. 

Proof. First, the new constructed method (2) is applied to the test equation 


(14), and then both sides of the obtained equation are collocated at points 


cq, t=1,...,m. Thus the algebraic relation is obtained in the form 
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P(tn + ch) = Svat (ci) Yn- ME wle) (c;)P (ta,7) ) + A*h? Y wile) (ci) P(tnj); 


(17) 


fori =1,...,m. Also, the approximated solution in the next step is obtained 
by 


r-1 
inti =>) oe ane +n v4 (1) Pltng) +R? >) xi(L)Pltn,3), (18) 
k=0 j=l j=l 


forn =r—1,r,...,N-—1. By the notations introduced in (15), the relations 


(17)-(18) can be written in the matrix form as 


P(tn + ch) = v(c)y? + zAP(tn + ch) + 27AP(ty + ch), 


1 
Yn+1 = y(1)y” + zu? P(ty +ch) +z 2wt P(tn + ch). a 
Thus by setting Q =I — zA — z?A, the first relation in (19) leads to 
Pltn + ch) = Q"*y(o)y®”, 
where substituting it in (19) for computing y,+1 leads to 
Yass = (ell) + (20? + 2?wT)Q7*9(c))¥$. (20) 


This relation is equivalent by the first row of (16) and the proof is complete. 


The coefficient matrix in (16) is shown by R(z) € C&T) *("+)) and it is 
called the stability matrix of the method. So the stability function is obtained 
by 

p(w, z) = det(wI — R(z)). (21) 


The stability properties of (2) with respect to the standard test equation (14) 
are dependent on the corresponding stability function p(w, z). 

Now, by multiplying the stability function (21) by its denominator, the 
stability polynomial of the method is obtained, which will be shown by the 
same symbol p(w, z). Therefore, the stability properties of the corresponding 


methods depend on the obtained polynomial p(w, z). The region of absolute 
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stability, A, is defined to be a region of the complex z-plane such that the 
roots of p(w, z) = 0 lie within the unit circle whenever z lies in the interior 
of the region. To obtain the region of absolute stability, we use the bound- 
ary locus method [15]. Inserting w = e’, the roots of stability polynomial 


describe the bound of stability region. 


5 Construction of A-stable SDMCMs 


In this section, we focus our attentions on two steps and three steps con- 
tinuous methods, (r = 2,3) with one and two collocation parameters. The 
polynomials of the methods are constructed by solving the linear system (8). 
Also, one can construct these polynomials by solving the system of order 
conditions (11) by p = 2m+r-—1. Then, by considering the stability matrix 
(16), collocation parameters in the interval [0,1] are determined in order to 


achieve A-stability. 


5.1 Construction of SDMCMs with r = 2 and m=1 


In this subsection, we first investigate two step new methods with r = 2 and 
one collocation parameter of order p = 2m +r—1 = 3. The coefficients of 


these methods are 


s? — 3cs? + 3¢?s + 3c? ++ 3c4+1 


yo(s) = 5 Se 44 ’ 
yi(s) = le cee 

010) 

ei s((2c + 1)s? — (3c? +1)s — 3c? 2c) 


2(3c? + 3c + 1) 


The stability polynomial of the method is written in the form 
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By performing an advanced search based on the boundary locus method, we 
find that there are some A-stable methods within this class of methods, for 
example, when the collocation parameter c lies in the interval [0.578, 1], the 
obtained method is the A-stable of order 3. For example, by choosing c = 1, 
an A-stable method of order 3 is obtained where the polynomials of methods 


are in the form 


s? — 387 +38+7 —s® + 3s? — 3s 

Yo(s) = ) yi(s) —  —_ oo - os 
7 7 

—s® + 3s? + 4s 35° — 2s? — 5s 

t(s) = =F) as) = 


5.2 Construction of SDMCMs with r = 2 and m = 2 


In this subsection, we describe the construction of the method with r = 2 
and two collocation parameters. The polynomials of method are of degree 
5 and can be obtained by solving the system of order condition (11) with 
p=2m+r—1=5. The stability polynomial of the method is obtained in 


the form 
7 * 
p(w, z) = So pi(zyw', 
t=1 


where p;(z) are polynomials of degree 2. By an extensive computer searching, 
we could not found A-stable methods for c1,c2 € [0,1], but in some cases, an 


extensive stability region (near to A-stability) can be observed. 


For example, choosing c = [5 1] leads to method of order 5 with un- 
bounded stability region, where its stability region is plotted in Figure 1, 
and its coefficients are determined by 

15 455,53 454, 65 
182” 182° "14° 182° ' 91°’ 
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p(s) = = oe me ms im * mi oh 
Oe, he a yy 
i= Rte Me he Bae 
gg a a oa 


Also, the stability polynomial of this method is 


p(w,z) = x45 (1824 — 15023 + 7272? — 21202 + 2912)u? 
t atts ((—722? — 7682 — 2944)w + 2? + 8z + 32). 


real(z) 


Figure 1: Stability region (Shaded region) of method with r= 2, m = 2 and ¢ = [5, 1]. 


5.3 Construction of SDMCMs with r= 3 and m=1 


We now consider 3-step SDMCMs with one collocation parameter of order 
p =A. Solving the system of order condition (11) corresponding to p = 4, 
we obtain the family of methods of order 4 depending on the parameter c. 


We apply the boundary locus method on the stability polynomial p(w, z), in 
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order to determine the values of the free parameters c € [0,1] achieving A- 
stability. We obtained when c belongs in [0.619, 1], the constructed methods 
are A-stable. Choosing, for example, for c = i. we obtain the three-step 


formula of uniform order p = 4 with coefficients given by 


5.4 Construction of SDMCMs with r = 3 and m = 2 


We now present the construction for SDMCMs of uniform convergence order 
6 with r = 3 and m = 2. First, by solving the order conditions (11), the 
polynomials of a method of degree 6 are determined depending on the pa- 
rameters cy and cg. The stability polynomial of the method is obtained in 


the form 
7 * 
p(w, z) = So pi(z)w', 
t=1 


where p;(z) are polynomials of degree 2. We apply the boundary locus 
method on the stability polynomial p(w, z), in order to determine the val- 
ues of the free parameters c1,c2 € [0,1] achieving A-stability. The range 
of (c1,¢2) in the domain [0,1] x [0,1], which leads to A-stable methods, is 
plotted in Figure 2. 


For example, by choosing c = [5 1], an A-stable method of order p = 6 is 
obtained, where coefficients are 


_ 13216s° — 21564s° — 33123s4 + 1012598? — 87639s + 32517s + 219758 


~o(s) 219758 : 
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146728 + 25536s° + 30786s* — 103768s° + 913088? — 34104s 


vis) = 219758 
) 1456s® — 3972s° + 2337s* + 2509s? — 3669s? + 1587s 
Ss = + I 
_ 219758 
ia) 77504s° + 20400s° — 305456s* + 707685? + 2538728? — 65248s 
s a % 
109879 
bo(s) 83384s® — 11604s° + 323186s* — 120143s* — 2118878? + 159662s 
s ad + J 
: 109879 
) 2744s® — 680s° — 1112884 + 8296s? + 7520s? — 8480s 
Ss = ? 
oe 9989 
) 33320s® + 12944s° — 116167s* + 338728? + 760258? — 536388 
X2\s) = ‘ 


219758 


0.8 


0.4 


0.2 


0 02 04 0.6 08 1 


Figure 2: Acceptable (ci, c2) pairs for A-stability of SDMCMs with r = 3, m = 2. 


6 Numerical eexperiments 


In this section, we present numerical results arising from the application 
of the SDMCMs in order to confirm the theoretical expectations, show the 
efficiency and accuracy of the new method, and validate the order of these 
methods in the integration of stiff systems. In what follows, we describe the 


details of the implemented methods: 


e Method 1: A-stable SDMCM of convergence order 3 with r = 2, m = 


1, and collocation parameter c = 1. 
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e Method 2: SDMCM of convergence order 5 with r = 2, m = 2, and 


collocation parameters c, = on c@=l1. 


e Method 3: A-stable SDMCM of convergence order 4 with r = 3, m= 


1, and collocation parameter c = 1. 


e Method 4: A-stable SDMCM of convergence order 6 with r = 3, m= 


2, and collocation parameters c = [4,1]. 


Also, we compared the results with the method in [6, 8], by the implementa- 
tion of the method. 


e Method 5: Second derivative two-step collocation method with m = 1, 


p=4m=4,c= 5 and gq = 1, introduced in [8]. 


e Method 6: Two-step almost collocation method with m = 2, p = 


Im 4, e= [2 


5) ap) introduced in [6]. 


Computational experiments are done by applying the methods to the follow- 
ing problems: 
P1. The nonlinear stiff ODE 

y(t) = —10004y, (t) + 10000y2(t)?, yi (0) = 1, 


yo(t) = yr(t) — yo(t)(1+ yo(t)®), — y2(0) = 1, 


with t € [0,1]. The exact solution is [y:(t) yo(t)|? = [e~*® e74]". 
P2. The Robertson chemical kinetics problem [13] 
y(t) = —0.04y1 (t) + 10*y2(t)ys(t), yi(0) = 1, 
Yo(t) = 0.04y; (t) — 10*y2(t)ya(t) — 3 x 10"ye(t)?, yo(0) = 0, 
Y3 = 3 x 10% yo(t)?, y3(0) = 0, 


with t € [0, 1000]. 


P3. The Pleiades problem [16] is a celestial mechanics problem of seven 
stars in the plane of coordinates x; and y; and masses m; = i for i = 
1,2,...,7. The problem consists of a system of 14 special second-order 


differential equations rewritten to a first-order form, thus providing a 
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system of ODEs of dimension 28. The formulation and data have been 


taken from [11]. The problem is in the form 
2" = f(z), 20) =20, 2'(0) = 2, 


with z € R'4 and 0 < t < 3. By considering z = (27, y7)’, 
z,y € R’, the function f : R'* > R"* is given by f(z) = f(z,y) = 
(f(x,y), f(x, y))”, where f, f@) : RM — R? are 


£P =o yay — 0) (ri, 


i#i 
2 3/2 - 
fp = So my(y; — y)/r3f , t=1,2,...,7, 
i#i 
with m; =i and 
Ti = (x; — 23)? + (yi = wr 


We convert this problem to a first-order form by defining w = z’, which 


leads to a system of 28 nonlinear differential equations of the form 


where (z7,w7)? € R*8 and 0 <t < 3. The initial values are 


3, 3, —i, —3, 2, —2, ay, 


xo ro = ( 
zo \ | Yo yo = (3, —3, 2,0, 0, —4,4)7, 
(:) ~ lat 1?) af = (0,0,0,0,0, 1.75, -1,5)*, 
yy yt = (0,0, 0, —1.25, 1,0,0)7. 


Table 1 presents the reference solution at the end of the integration 


interval, which is reported in [16]. 


In our numerical experiments, we apply Method 1— Method 4 to the given test 
problems with the fixed step sizes h = T/N for several integer values of N. 
The global error of the methods at the endpoint of the interval of integration 
is listed by Err. Also, to verify theoretical results on the order of accuracy, 
we compute a numerical estimate to the order of accuracy of the methods 


by the formula log, (||Errn(x)||/||Errn2(x)||), where Errp(x) and Err; 2(2) 
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Table 1: Reference solution at the end of the integration interval 


x1 =0.3706139143970502 


y1 = -3.943437585517392 


£2 = 3.237284092057233 


y2 = -3.271380973972550 


%3 = -3.222559032418324 


y3 =5.225081843456543 


x4 = 0.6597091455775310 


ya =-2.590612434977470 


5 = 0.3425581707156584 


ys = 1.198213693392275 


tg=1.562172101400631 


Ye = -0.2429682344935824 


£7 = -0.7003092922212495 


Y7 = 1.091449240428980 


x’, = 3.417003806314313 


y, = —3.741244961234010 


x = 1.354584501625501 


Yo = 0.3773459685750630 


x’; = —2.590065597810775 


Y3 = 0.9386858869551073 


x, = 2.025053734714242 


y4 = 0.3667922227200571 


x, = —1.155815100160448 


ys = —0.3474046353808490 


x = —0.8072988170223021 


YG = 2.344915448180937 


xt = 0.5952396354208710 


ys = —1.947020434263292 


386 


are the errors at the endpoint of the interval of integration corresponding to 


the step sizes h and h/2, respectively. we observed that the expected order 


was achieved. Also, by comparing results by methods introduced in [6, 8], 


one can see that the results are compatible with Method 5 and the solution 


by this new method is more precise than the results of Method 4 while the 


computational cost for the new method is less. 


In problems P2, the reference solutions reported in [10] are used for com- 


puting the global error. The obtained results are reported in Table 2. 
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Table 2: The results of the methods for problem P1. 


N 4 8 16 32 64 

Method1 = Err | 2.22x10-4 3.40x 10-5 4.64x 10-© =6.04x 10-7 ~—_7.71 x 10-8 
p(h) 2.71 2.87 2.94 2.97 

Method 2. Err | 2.12x 10-7 8.38x10-® 2.93x10-! 966x107! 3.10x 10-8 
p(h) 4.66 4.84 4,92 4.96 

Method 3. Err | 1.95x 1075 1.92x10-® 1.39x1077 9.301079 5.99 x 1071 
p(h) 3.35 3.78 3.91 3.96 

Method 4 Err | 3.96 x 1079 1.01x107'® 1.93x 1071? 3.33x 1074 5.46 x 10716 
p(h) 5.31 5.70 5.86 5.93 

Method 5 Err | 1.17x 10-3 6.14x 10-5 = 2.58x10-® 1.01x 10-7 5.96 x 10-9 
p(h) 4.26 4.57 4.67 4.59 

Method 6 Err | 1.51x 1073 =1.74x 10-4 = 1.49x 10-® =: 1.08 x 10° ~—- 7.08 x 10-8 
p(h) 3.11 3.54 3.79 3.93 
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Table 3: The results of SDMCMs for problem P2 


N 500 750 1000 50 1500 

Method 1 Err | 4.41x 107° 1.69x 10-® 8.22x107-® 462x107® 2.85 x 1077 
p(h 2.37 2.50 58 2.69 

Method 2. Err | 3.41107" 7.41x 10-8 2.17x107-8 8.07x10-9 4.46 x 10-9 
p(h 3.76 4.27 4.42 4.64 

Method 3 Err | 1.99x10-® 8.24x10-7 399x107? 216x10-* 7.81 x 1078 
p(h 3.07 3.25 3.36 3.54 

Method 4 Err | 7.861078 1.67x 107-8 5.261079 1.5 0-8 5.42 x 10719 


4.84 


4.97 


5.37 5.90 


Table 4: The results of SDMCMs for problem P3 


h 4 x 10-8 4 Xx 1078 x x 1073 ig X 1073 
Method1 = Err | 199x107! 2.49x 10-2 313x107? 3.90 x 1074 
p(h) 2.99 2.99 3.00 
Method 2 = Err | 7.32x10-* 2.31x 1075 7.21x107-7 2.47 x 10-8 
p(h) 4.98 4.99 5.00 
Method 3. Err | 5.39x 107-2 2.52x 107% 1.24x10-* 6.67 x 10-4 
p(h) 4.41 4.34 4,22 
Method 4 = Err | 2.15x 107-5 3.71x 107-7 6.02x10-9 9.46 x 1071! 
p(h) 5.86 5.94 5.99 


7 Conclusion 


We have introduced a new family of SDMCMs. 


The constructed r-step 


method with m collocation parameters has a uniform order 2m +r — 1. 


Examples of SDMCMs up to order 6 with the property of A-stability were 


constructed. The accuracy and efficiency of constructed methods were ver- 


ified by solving some stiff problems. SDMCMs are efficient in solving stiff 
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problems, which are confirmed by numerical experiments. Future work will 
be in the investigation of G-simplecticity of SDMCMs [12, 17]. 
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